﻿using System;
using System.Collections.Generic;

namespace LxCFD
{
    public partial class VelocityVector : MassFluxBase
    {
        public List<NoSlipWall> NoSlipWalls { get; set; } = new List<NoSlipWall>();
        public List<Symmetry> Symmetries { get; set; } = new List<Symmetry>();
        public List<VelocityInlet> VelocityInlets { get; set; } = new List<VelocityInlet>();
        public List<PressureInOut> PressureInOuts { get; set; } = new List<PressureInOut>();
        public Value1Base<double> UVW_Gamma { get; set; }
        public ValueBase<Vector3d> UVW_Source;
        public Value1Base<Vector3d> UVW_IniValue;
        public double UVW_RelaxFactor = 1;
        public Value1Base<double> DynamicViscosity;

        public NSEquation U;
        public NSEquation V;
        public NSEquation W;
        protected double[] U_x;
        protected double[] V_x;
        protected double[] W_x;

        protected double[] U_bFaceCommon;
        protected double[] V_bFaceCommon;
        protected double[] W_bFaceCommon;

        public Value1Base<double> Density;

        public Vector3d GetVelocity(int cellIndex)
        {
            return new Vector3d(U_x[cellIndex], V_x[cellIndex], W_x[cellIndex]);
        }
        public override void CheckCreate(FvZone zone)
        {
            base.CheckCreate(zone);

            this.U_bFaceCommon = zone.FaceCommonArray;
            LxTool.CheckCreateArray(ref V_bFaceCommon, Faces.Length);
            LxTool.CheckCreateArray(ref W_bFaceCommon, Faces.Length);

            if (U == null)
            {
                U = new NSEquation();
                U.CheckCreate(zone);
            }
            this.U_x = U.x;
            U.b_FaceCommon = this.U_bFaceCommon;
            if (V == null)
            {
                V = new NSEquation();
                V.CheckCreate(zone);
            }
            this.V_x = V.x;
            V.b_FaceCommon = this.V_bFaceCommon;
            if (W == null)
            {
                W = new NSEquation();
                W.CheckCreate(zone);
            }
            this.W_x = W.x;
            W.b_FaceCommon = this.W_bFaceCommon;
        }
        public void SetZone_CommonTerm(Value1Base<double> UVW_gamma, Value1Base<double> dynamicViscosity, ValueBase<Vector3d> UVW_source, Value1Base<Vector3d> UVW_iniValue, double UVW_relaxFactor)
        {
            NoSlipWalls.Clear();
            Symmetries.Clear();
            VelocityInlets.Clear();
            PressureInOuts.Clear();

            this.UVW_Gamma = UVW_gamma;
            this.DynamicViscosity = dynamicViscosity;
            this.UVW_Source = UVW_source;
            this.UVW_IniValue = UVW_iniValue;
            this.UVW_RelaxFactor = UVW_relaxFactor;

            U.SetZone_CommonTerm(this.UVW_Gamma, Value1Const.Double0, null, this.UVW_RelaxFactor);
            V.SetZone_CommonTerm(this.UVW_Gamma, Value1Const.Double0, null, this.UVW_RelaxFactor);
            W.SetZone_CommonTerm(this.UVW_Gamma, Value1Const.Double0, null, this.UVW_RelaxFactor);

        }
        public void SetZone_ConvTerm(Value1Base<double> density, ConvectionSchemes convectionSchemes)
        {
            this.Density = density;
            Value1Const.Double conv1 = new Value1Const.Double(1);
            U.SetZone_ConvTerm(conv1, convectionSchemes, this);
            V.SetZone_ConvTerm(conv1, convectionSchemes, this);
            W.SetZone_ConvTerm(conv1, convectionSchemes, this);
        }
        public void SetBoundary_NoSlipStaticWall(FvThread thread)
        {
            SetBoundary_NoSlipMovingWall(Value1Const.Vector0, thread);
        }
        public void SetBoundary_NoSlipMovingWall(Value1Base<Vector3d> boundaryVelocity, FvThread thread)
        {
            NoSlipWall noSlipWall = new NoSlipWall();
            noSlipWall.SetHead(thread);
            noSlipWall.SetParameter(boundaryVelocity);
            NoSlipWalls.Add(noSlipWall);

            U.SetBoundary_symmetry(thread);
            V.SetBoundary_symmetry(thread);
            W.SetBoundary_symmetry(thread);
        }
        public void SetBoundary_SlipWall(FvThread thread)
        {
            Symmetry boundarySlipWall = new Symmetry();
            boundarySlipWall.SetHead(thread);
            Symmetries.Add(boundarySlipWall);

            U.SetBoundary_symmetry(thread);
            V.SetBoundary_symmetry(thread);
            W.SetBoundary_symmetry(thread);
        }
        public void SetBound_VelocityInlet(Value1Base<double> value, FvThread thread)//方向垂直壁面指向域内
        {
            SetBound_VelocityInlet(new Value3dBoundInVel(value, Zone), thread);
        }
        public void SetBound_VelocityInlet(Value1Base<Vector3d> value3d, FvThread thread)
        {
            VelocityInlet velinlet = new VelocityInlet();
            velinlet.SetHead(thread);
            velinlet.SetParameter(value3d);
            VelocityInlets.Add(velinlet);

            U.SetBoundary_flow(null, thread);
            V.SetBoundary_flow(null, thread);
            W.SetBoundary_flow(null, thread);
        }
        public void SetBound_PressureInOut(Value1Base<double> pressure, FvThread thread)//方向垂直壁面指向域内
        {
            Value3dBoundNormalDir dir = new Value3dBoundNormalDir(Zone);
            SetBound_PressureInOut(pressure, dir, thread);
        }
        public void SetBound_PressureInOut(Value1Base<double> pressure, Value1Base<Vector3d> direction, FvThread thread)
        {
            PressureInOut pressureinout = new PressureInOut();
            pressureinout.SetHead(thread);
            pressureinout.SetParameter(pressure, direction);
            PressureInOuts.Add(pressureinout);

            U.SetBoundary_flow(null, thread);
            V.SetBoundary_flow(null, thread);
            W.SetBoundary_flow(null, thread);
        }
        #region 外部调用
        public override void Initialize()
        {
            MultiThreads.RunCalTasks(initialize);
        }
        public void initialize(int taskIndex)
        {
            setUVW_x(UVW_IniValue, InnerIndex, taskIndex);
            SetUVW_x(taskIndex);
            foreach (var b in Symmetries)
            {
                set_x_SlipWall(b.ThreadIndex, taskIndex);
            }

            foreach (var b in PressureInOuts)
            {
                setUVW_x(UVW_IniValue, b.ThreadIndex, taskIndex);
            }
            //重置mass相关
            setFaceFlux_0(Zone.TotalIndex, taskIndex);
            setMassImbanlance_0(Zone.TotalIndex, taskIndex);
        }
        public void SetUVW_x(int taskIndex)
        {
            foreach (var b in NoSlipWalls)
            {
                if (b.SetVelState.GetIsNeedToCal(taskIndex))
                {
                    if (b.boundaryVelocity.IsEqual0())//(b.boundaryVelocity is Value3dConst vel) && vel.vector3d.Is0()
                    {
                        setUVW_x(b.boundaryVelocity, b.ThreadIndex, taskIndex);
                    }
                    else
                    {
                        set_x_MovingWall(b.boundaryVelocity, b.ThreadIndex, taskIndex);
                    }
                }
            }
            foreach (var b in VelocityInlets)
            {
                if (b.SetVelState.GetIsNeedToCal(taskIndex))
                {
                    setUVW_x(b.Velocity, b.ThreadIndex, taskIndex);
                }
            }
        }
        public void AddViscositySource_LoopFace(int taskIndex)//由于每个线程处理的index是固定的，不需要考虑同步的问题
        {
            //由于粘性导致的Source
            addFaceViscositySource_inner(taskIndex);
            foreach (var b in NoSlipWalls)
            {
                addUVW_Ab_WallShearStress(b.ThreadIndex, taskIndex);//不采用fixValue的形式，可以实现层流与湍流的通用性，按理说还得添加垂直壁面方向的力，没有加
                addUVW_Ab_NormalStress(DynamicViscosity, b.ThreadIndex, taskIndex);
            }
            foreach (var b in Symmetries)
            {
                addUVW_Ab_NormalStress(UVW_Gamma, b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityInlets)
            {
                addFaceViscositySource_boundOwner(b.ThreadIndex, taskIndex);
            }
            foreach (var b in PressureInOuts)
            {
                addFaceViscositySource_boundOwner(b.ThreadIndex, taskIndex);
            }
        }
        public void AddSource_LoopCell(int taskIndex)
        {
            //附加的Source,Scalar中的Source
            Add_bSource(UVW_Source, InnerIndex, taskIndex);
        }
        public void UpdateUVWBound(int taskIndex)
        {
            foreach (var b in Symmetries)
            {
                set_x_SlipWall(b.ThreadIndex, taskIndex);
            }
            foreach (var b in PressureInOuts)
            {
                set_x_pressureBound(b.Direction, b.ThreadIndex, taskIndex);
            }
#if DEBUG//仅调试用
            for (int i1 = InnerIndex.FaceIDs.End + 1; i1 <= TotalIndex.FaceIDs.End; i1++)
            {
                ref FvFace face = ref Faces[i1];
                MassImbalance[face.NSideCell] = MassImbalance[face.OSideCell];
            }
#endif
        }

        #endregion
        #region 内部UVW相关方法
        protected void setUVW_x(Value1Base<Vector3d> UVWValue, IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                Vector3d v = UVWValue.GetCellValue(i1);
                U_x[i1] = v.X;
                V_x[i1] = v.Y;
                W_x[i1] = v.Z;
            }
        }
        protected void set_x_MovingWall(Value1Base<Vector3d> UVWValue, IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d v = UVWValue.GetCellValue(face.NSideCell);
                v = Vector3d.VectorToPlane(v, face.Normal);
                U_x[face.NSideCell] = v.X;
                V_x[face.NSideCell] = v.Y;
                W_x[face.NSideCell] = v.Z;
            }
        }
        protected void set_x_SlipWall(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d v = new Vector3d(U_x[face.OSideCell], V_x[face.OSideCell], W_x[face.OSideCell]);
                v = Vector3d.VectorToPlane(v, face.Normal);
                U_x[face.NSideCell] = v.X;
                V_x[face.NSideCell] = v.Y;
                W_x[face.NSideCell] = v.Z;
            }
        }
        public void set_x_pressureBound(Value1Base<Vector3d> direction, IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                //要区分流出还是流入
                double massflux = MassFlux[i1];
                if (massflux > 0)////由明确的流出时，确定为流出边界
                {
                    U.CalBoundValue_grad0(ref face);
                    V.CalBoundValue_grad0(ref face);
                    W.CalBoundValue_grad0(ref face);
                }
                else
                {
                    Vector3d dir = direction.GetFaceValue(i1);
                    double des = Density.GetFaceValue(i1);
                    double evsb = Vector3d.DotProduct(dir, face.Normal);
                    if (evsb != 0)
                    {
                        double massb = massflux;
                        Vector3d v = massb / (des * evsb) * dir;
                        U_x[face.NSideCell] = v.X;
                        V_x[face.NSideCell] = v.Y;
                        W_x[face.NSideCell] = v.Z;
                    }
                    else
                    {
                        LxConsole.Error("压力边界所指定的速度方向不指向区域内");
                    }
                }
            }
        }
        protected void Add_bSource(ValueBase<Vector3d> S_Value3d, IndexNumber indexNumber, int taskIndex)
        {
            double[] U_ADiag = U.ADiag; double[] U_b = U.b;
            double[] V_ADiag = V.ADiag; double[] V_b = V.b;
            double[] W_ADiag = W.ADiag; double[] W_b = W.b;
            if (S_Value3d is Value1Base<Vector3d> s3d1)
            {
                if (s3d1.IsEqual0()) return;//s3d1 is Value3dConst vc && vc.vector3d.Is0()
                indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    Vector3d v = s3d1.GetCellValue(i1);
                    v.MultiplyToSelf(Cells[i1].Volumn);
                    U_b[i1] += v.X;
                    V_b[i1] += v.Y;
                    W_b[i1] += v.Z;
                }
            }
            else if (S_Value3d is Value2Base<Vector3d> s3d2)
            {
                indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    double volumn = Cells[i1].Volumn;
                    s3d2.GetCellValue(out Vector3d va, out Vector3d vb, i1);
                    va.MultiplyToSelf(volumn);
                    vb.MultiplyToSelf(volumn);
                    U_b[i1] += vb.X;
                    V_b[i1] += vb.Y;
                    W_b[i1] += vb.Z;
                    U_ADiag[i1] -= va.X;
                    V_ADiag[i1] -= va.Y;
                    W_ADiag[i1] -= va.Z;
                }
            }
            else
            {
                LxConsole.Error(string.Format($"区域{Zone.name} 速度矢量场无法识别源项"));
            }
        }
        protected void addFaceViscositySource_inner(int taskIndex)
        {
            Vector3d[] U_grad = U.Grad; double[] U_bFaceCommon = U.b_FaceCommon;
            Vector3d[] V_grad = V.Grad; double[] V_bFaceCommon = V.b_FaceCommon;
            Vector3d[] W_grad = W.Grad; double[] W_bFaceCommon = W.b_FaceCommon;
            InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d dU = face.CalCFCorrectedGradient_inner(U_grad, U_x);
                Vector3d dV = face.CalCFCorrectedGradient_inner(V_grad, V_x);
                Vector3d dW = face.CalCFCorrectedGradient_inner(W_grad, W_x);
                double gamma = UVW_Gamma.GetFaceValue(i1);
                U_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.X, dV.X, dW.X), face.Normal);
                V_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.Y, dV.Y, dW.Y), face.Normal);
                W_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.Z, dV.Z, dW.Z), face.Normal);
            }
        }
        protected void addFaceViscositySource_boundOwner(IndexNumber indexNumber, int taskIndex)
        {
            Vector3d[] U_grad = U.Grad; double[] U_bFaceCommon = U.b_FaceCommon;
            Vector3d[] V_grad = V.Grad; double[] V_bFaceCommon = V.b_FaceCommon;
            Vector3d[] W_grad = W.Grad; double[] W_bFaceCommon = W.b_FaceCommon;
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d dU = face.CalCFCorrectedGradient_bound(U_grad, U_x);
                Vector3d dV = face.CalCFCorrectedGradient_bound(V_grad, V_x);
                Vector3d dW = face.CalCFCorrectedGradient_bound(W_grad, W_x);
                double gamma = UVW_Gamma.GetFaceValue(i1);
                U_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.X, dV.X, dW.X), face.Normal);
                V_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.Y, dV.Y, dW.Y), face.Normal);
                W_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.Z, dV.Z, dW.Z), face.Normal);
            }
        }
        protected void addUVW_Ab_WallShearStress(IndexNumber indexNumber, int taskIndex)
        {
            double[] U_ALUN = U.ALUN; double[] U_bFaceCommon = U.b_FaceCommon;
            double[] V_ALUN = V.ALUN; double[] V_bFaceCommon = V.b_FaceCommon;
            double[] W_ALUN = W.ALUN; double[] W_bFaceCommon = W.b_FaceCommon;
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d ov = GetVelocity(face.OSideCell);
                Vector3d nv = GetVelocity(face.NSideCell);
                double XX = Math.Pow(face.Normal.X, 2);
                double YY = Math.Pow(face.Normal.Y, 2);
                double ZZ = Math.Pow(face.Normal.Z, 2);
                double XY = face.Normal.X * face.Normal.Y;
                double YZ = face.Normal.Y * face.Normal.Z;
                double XZ = face.Normal.X * face.Normal.Z;
                double commonterm = UVW_Gamma.GetFaceValue(i1) / Vector3d.DotProduct(face.OToN, face.Normal);
                double temp = YY + ZZ;
                U_ALUN[i1] -= commonterm * (temp);
                U_bFaceCommon[i1] += commonterm * (
                     nv.X * (temp)
                    + (ov.Y - nv.Y) * XY
                    + (ov.Z - nv.Z) * XZ
                    );
                temp = XX + ZZ;
                V_ALUN[i1] -= commonterm * (temp);
                V_bFaceCommon[i1] += commonterm * (
                      (ov.X - nv.X) * XY
                    + nv.Y * (temp)
                    + (ov.Z - nv.Z) * YZ
                    );
                temp = XX + YY;
                W_ALUN[i1] -= commonterm * (temp);
                W_bFaceCommon[i1] += commonterm * (
                      (ov.X - nv.X) * XZ
                    + (ov.Y - nv.Y) * YZ
                    + nv.Z * (temp)
                    );
            }
        }
        protected void addUVW_Ab_NormalStress(Value1Base<double> gamma, IndexNumber indexNumber, int taskIndex)
        {
            double[] U_ALUN = U.ALUN; double[] U_bFaceCommon = U.b_FaceCommon;
            double[] V_ALUN = V.ALUN; double[] V_bFaceCommon = V.b_FaceCommon;
            double[] W_ALUN = W.ALUN; double[] W_bFaceCommon = W.b_FaceCommon;
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d ov = GetVelocity(face.OSideCell);
                double common = 2 * gamma.GetFaceValue(i1) / Vector3d.DotProduct(face.OToN, face.Normal);
                U_ALUN[i1] -= common * Math.Pow(face.Normal.X, 2);
                U_bFaceCommon[i1] -= common * (ov.Y * face.Normal.Y + ov.Z * face.Normal.Z) * face.Normal.X;
                V_ALUN[i1] -= common * Math.Pow(face.Normal.Y, 2);
                V_bFaceCommon[i1] -= common * (ov.X * face.Normal.X + ov.Z * face.Normal.Z) * face.Normal.Y;
                W_ALUN[i1] -= common * Math.Pow(face.Normal.Z, 2);
                W_bFaceCommon[i1] -= common * (ov.X * face.Normal.X + ov.Y * face.Normal.Y) * face.Normal.Z;
            }
        }
        #endregion
        #region boundaryTypes
        public class Wall : BoundaryHeadBase
        {
            public ValueStateLog FluxState { get; set; } = new ValueStateLog(true);
        }
        public class NoSlipWall : Wall
        {
            public Value1Base<Vector3d> boundaryVelocity;
            public ValueStateLog SetVelState { get; set; }
            public void SetParameter(Value1Base<Vector3d> boundaryVelocity)
            {
                this.boundaryVelocity = boundaryVelocity;
                SetVelState = new ValueStateLog(boundaryVelocity.IsConstant());
            }
        }
        public class Symmetry : Wall
        {
        }
        public class VelocityInlet : BoundaryHeadBase
        {
            public Value1Base<Vector3d> Velocity { get; private set; }
            public ValueStateLog SetVelState { get; set; }
            public void SetParameter(Value1Base<Vector3d> velocity)
            {
                this.Velocity = velocity;
                SetVelState = new ValueStateLog(velocity.IsConstant());
            }
        }
        public class PressureInOut : BoundaryHeadBase
        {
            public Value1Base<double> Pressure { get; private set; }
            public Value1Base<Vector3d> Direction { get; private set; }//不需要为单位向量
            public ValueStateLog SetPressureState;
            public void SetParameter(Value1Base<double> Pressure, Value1Base<Vector3d> Direction)
            {
                this.Pressure = Pressure;
                this.Direction = Direction;
                SetPressureState = new ValueStateLog(Pressure.IsConstant());
            }
        }
        public class Value3dBoundInVel : Value1Base<Vector3d>
        {
            public Value3dBoundInVel(Value1Base<double> value, FvZoneBase zone)
            {
                this.value = value;
                this.Faces = zone.Faces;
                this.Cells = zone.Cells;
            }
            FvFace[] Faces;
            FvCell[] Cells;
            Value1Base<double> value;
            public override Vector3d GetCellValue(int cellIndex)
            {
                ref FvFace face = ref Faces[Cells[cellIndex].NFaces[0]];
                return (-value.GetCellValue(cellIndex) / face.Area) * face.Normal;
            }
        }
        public class Value3dBoundNormalDir : Value1Base<Vector3d>
        {
            public Value3dBoundNormalDir(FvZoneBase zone)
            {
                this.Faces = zone.Faces;
            }
            FvFace[] Faces;
            public override Vector3d GetFaceValue(int faceIndex)
            {
                return Faces[faceIndex].Normal;
            }
        }
        #endregion
    }
    public partial class VelocityVector : ICalResidual
    {
        protected double Max_Rc;
        protected double Max_acPhic;
        public void CalResidual(int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double maxFaceFlux = 0, maxImbalance = 0;
            for (; i1 < endp1; i1++)
            {
                if (MassImbalance[i1] > maxImbalance) maxImbalance = MassImbalance[i1];
            }
            TotalIndex.FaceIDs.GetSeg(out i1, out endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                if (MassFlux[i1] > maxFaceFlux) maxFaceFlux = MassFlux[i1];
            }
            ThreadsMath.KeepMax(ref this.Max_Rc, maxImbalance);
            ThreadsMath.KeepMax(ref this.Max_acPhic, maxFaceFlux);
        }
        public void GetResidual_Reset(out double Rc, out double acPhic)
        {
            Rc = this.Max_Rc;
            acPhic = this.Max_acPhic;
            this.Max_Rc = 0;
            this.Max_acPhic = 0;
        }
    }

    public partial class VelocityVector//界面速度插值及massFlux计算
    {
        protected RhieChowInterpolation VelocityInterpolation;
        public override void CalFaceMassFlux(int taskIndex)
        {
            calFaceMassFlux_inner(taskIndex);
            foreach (var b in NoSlipWalls)
            {
                setWallFlux0(b, taskIndex);
            }
            foreach (var b in Symmetries)
            {
                setWallFlux0(b, taskIndex);
            }
            foreach (var b in VelocityInlets)
            {
                calFaceMassFlux_VelBound(b.ThreadIndex, taskIndex);
            }
            foreach (var b in PressureInOuts)
            {
                calFaceMassFlux_pressureBound(b.ThreadIndex, taskIndex);
            }
        }
        private void setWallFlux0(Wall b, int taskIndex)
        {
            if (b.FluxState.GetIsNeedToCal(taskIndex))
            {
                setFaceFlux_0(b.ThreadIndex, taskIndex);
            }
        }
        protected void calFaceMassFlux_VelBound(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d vf = GetVelocity(face.NSideCell);//不进行修正，因为压力本身是根据内部网格计算出来的，不会出现棋盘型问题，且这样保证流入的质量一定等于速度乘以面积
                MassFlux[i1] = Density.GetFaceValue(i1) * Vector3d.DotProduct(vf, face.Normal);
            }
        }
        public void calFaceMassFlux_inner(int taskIndex)
        {
            InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d v_O = GetVelocity(face.OSideCell);
                Vector3d v_N = GetVelocity(face.NSideCell);
                Vector3d vf = (1 - face.gOwner) * v_O + face.gOwner * v_N + VelocityInterpolation.CalVelocityCorrection_Inner(i1, ref face);
                double des = Density.GetFaceValue(i1);
                MassFlux[i1] = des * Vector3d.DotProduct(vf, face.Normal);
                VelocityDirectionCorrection[i1] = calVelocityDirectionCorrecion(ref v_O, ref v_N, ref vf, ref face, des);
            }
        }
        public void calFaceMassFlux_pressureBound(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d v_O = GetVelocity(face.OSideCell);
                //Vector3d v_N = getVelocity(face.NSideCell);
                Vector3d vf = v_O + VelocityInterpolation.CalVelocityCorrection_Bound(i1, ref face);
                double des = Density.GetFaceValue(i1);
                MassFlux[i1] = des * Vector3d.DotProduct(vf, face.Normal);
                //VelDirCor应该不用修正，因为用不到
            }
        }
        protected float calVelocityDirectionCorrecion(ref Vector3d v_O, ref Vector3d v_N, ref Vector3d vf, ref FvFace face, double density)
        {
            //0为没有修正
            double vo2 = v_O.CalLengthSquared();
            double vn2 = v_N.CalLengthSquared();
            if (vo2 > 1e-100 && vn2 > 1e-100)
            {
                double cosa = Vector3d.DotProduct(v_O, v_N) / Math.Sqrt(vo2 * vn2);
                return (float)Math.Abs((1 - cosa) * density * Vector3d.DotProduct(vf, face.OToN));//* 0.5
            }
            else
            {
                return 0;// (float)(density.getFaceValue(faceIndex) * Vector3d.DotProduct(vf, face.OwnerToNeighbor));
            }
        }
    }
    public partial class VelocityVector  //瞬态处理
    {
        public TransientVelocity TransVelData;
        public void SetZone_TransTerm(bool isEnable, Value1Base<double> Trans, Value1Base<double> density)
        {
            U.SetZone_TransTerm(isEnable, Trans, density);
            V.SetZone_TransTerm(isEnable, Trans, density);
            W.SetZone_TransTerm(isEnable, Trans, density);
            if (isEnable)
            {
                if (TransVelData == null)
                {
                    TransVelData = new TransientVelocity();
                    TransVelData.Create(this);
                }
                TransVelData.SetTrans(density);
            }
            else
            {
                TransVelData = null;
            }
        }
        public void SetZone_RhieChow_AfterSetTrans(PressureCorrectionBasic pressureCorInSIMPLE, PressureSIMPLE pressureInSIMPLE)
        {
            if (TransVelData == null)
            {
                LxTool.CheckCreateInstance<RhieChowInterpolation, RhieChowInterpolation>(ref VelocityInterpolation);
                VelocityInterpolation.SetParameter(pressureCorInSIMPLE, pressureInSIMPLE);
            }
            else
            {
                TransVelData.PressureCor = pressureCorInSIMPLE;
                TransVelData.Pressure = pressureInSIMPLE;
                TransVelData.CheckCreateTransRhieChow();
            }
        }
        public class TransientVelocity : TransientDataBase
        {
            public Vector3d[] RhoVel_o;
            public Vector3d[] RhoVel_oo;
            public Vector3d[] Vel_o;
            public PressureCorrectionBasic PressureCor;
            public PressureSIMPLE Pressure;

            protected VelocityVector VelocityVector;
            protected Value1Base<double> Density;

            public void Create(VelocityVector velocityVector)
            {
                this.VelocityVector = velocityVector;
                LxTool.CheckCreateArray(ref RhoVel_o, velocityVector.Faces.Length);
                LxTool.CheckCreateArray(ref RhoVel_oo, velocityVector.Faces.Length);
                LxTool.CheckCreateArray(ref Vel_o, velocityVector.Faces.Length);
                //TransientDataState = TransientDataStates.None;
            }
            public void SetTrans(Value1Base<double> density)
            {
                this.Density = density;
            }
            public void CheckCreateTransRhieChow()
            {
                if (IsSecondOrder && TransientDataState == TransientDataStates.Both_o_oo)
                {
                    LxTool.CheckCreateInstance<RhieChow_2ndOrderTrans, RhieChowInterpolation>(ref VelocityVector.VelocityInterpolation);
                    var inter = VelocityVector.VelocityInterpolation as RhieChow_2ndOrderTrans;
                    inter.RhoVel_o = this.RhoVel_o;
                    inter.RhoVel_oo = this.RhoVel_oo;
                    inter.Vel_o = this.Vel_o;
                    inter.SetParameter(PressureCor, Pressure);
                }
                else
                {
                    LxTool.CheckCreateInstance<RhieChow_1stOrderTrans, RhieChowInterpolation>(ref VelocityVector.VelocityInterpolation);
                    var inter = VelocityVector.VelocityInterpolation as RhieChow_1stOrderTrans;
                    inter.RhoVel_o = this.RhoVel_o;
                    inter.Vel_o = this.Vel_o;
                    inter.SetParameter(PressureCor, Pressure);
                }
            }
            public override void Initialize()
            {
                MultiThreads.RunCalTasks(save0);
            }
            public override void PrepareToSave()
            {
                LxTool.Exchange(ref RhoVel_o, ref RhoVel_oo);
                switch (TransientDataState)
                {
                    case TransientDataStates.None: TransientDataState = TransientDataStates.Only_o; break;
                    case TransientDataStates.Only_o: TransientDataState = TransientDataStates.Both_o_oo; break;
                    default: break;
                }
                CheckCreateTransRhieChow();
            }
            public override void SaveOldValue(int taskIndex)
            {
                saveOldValue(VelocityVector.InnerIndex, taskIndex);
                foreach (var b in VelocityVector.PressureInOuts)
                {
                    saveOldValue(b.ThreadIndex, taskIndex);
                }
            }
            protected void saveOldValue(IndexNumber indexNumber, int taskIndex)
            {
                indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    RhoVel_o[i1] = Density.GetFaceValue(i1) * Vel_o[i1];
                }
            }
            protected void save0(int taskIndex)
            {
                VelocityVector.TotalIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    RhoVel_o[i1] = new Vector3d();
                }
            }
        }
    }
}

//U_bFaceCommon[i1] += commonterm * (
//      U_x[face.NSideCell] * (Math.Pow(face.Area, 2) - Math.Pow(face.Normal.X, 2))
//    + (V_x[face.OSideCell] - V_x[face.NSideCell]) * face.Normal.Y * face.Normal.X
//    + (W_x[face.OSideCell] - W_x[face.NSideCell]) * face.Normal.Z * face.Normal.X
//    );
//V_ALUN[i1] -= commonterm * (Math.Pow(face.Area, 2) - Math.Pow(face.Normal.Y, 2));
//V_bFaceCommon[i1] += commonterm * (
//      (U_x[face.OSideCell] - U_x[face.NSideCell]) * face.Normal.X * face.Normal.Y
//    + V_x[face.NSideCell] * (Math.Pow(face.Area, 2) - Math.Pow(face.Normal.Y, 2))
//    + (W_x[face.OSideCell] - W_x[face.NSideCell]) * face.Normal.Z * face.Normal.Y
//    );
//W_ALUN[i1] -= commonterm * (Math.Pow(face.Area, 2) - Math.Pow(face.Normal.Z, 2));
//W_bFaceCommon[i1] += commonterm * (
//      (U_x[face.OSideCell] - U_x[face.NSideCell]) * face.Normal.X * face.Normal.Z
//    + (V_x[face.OSideCell] - V_x[face.NSideCell]) * face.Normal.Y * face.Normal.Z
//    + W_x[face.NSideCell] * (Math.Pow(face.Area, 2) - Math.Pow(face.Normal.Z, 2))
//    );
//U_ALUN[i1] -= common * Math.Pow(face.Normal.X, 2);
//U_bFaceCommon[i1] -= common * (V_x[face.OSideCell] * face.Normal.Y + W_x[face.OSideCell] * face.Normal.Z) * face.Normal.X;
//V_ALUN[i1] -= common * Math.Pow(face.Normal.Y, 2);
//V_bFaceCommon[i1] -= common * (U_x[face.OSideCell] * face.Normal.X + W_x[face.OSideCell] * face.Normal.Z) * face.Normal.Y;
//W_ALUN[i1] -= common * Math.Pow(face.Normal.Z, 2);
//W_bFaceCommon[i1] -= common * (U_x[face.OSideCell] * face.Normal.X + V_x[face.OSideCell] * face.Normal.Y) * face.Normal.Z;

//public void RecordOldValue(int taskIndex)
//{
//    transientVelRecorder.RecordOldValue(taskIndex);
//}
//public void PrepareRecorder()
//{
//    transientVelRecorder.PrepareRecorder();
//    transientVelRecorder.CheckCreateTransRhieChow(ref velocityInterpolation);
//}
//public bool IsInitialized()
//{
//    return transientDataState != TransientDataStates.未记录;
//}
//public void SetInitialized()
//{
//    transientDataState = TransientDataStates.记录o;
//}
//foreach (var b in NoSlipWalls)
//{
//    if (b.setVelState.GetIsNeedToCal(taskIndex))
//    {
//        if (b.boundaryVelocity .isEqual0())//(b.boundaryVelocity is Value3dConst vel) && vel.vector3d.Is0()
//        {
//            setUVW_x(b.boundaryVelocity, b.indexNumber, taskIndex);
//        }
//        else
//        {
//            set_x_MovingWall(b.boundaryVelocity, b.indexNumber, taskIndex);
//        }
//    }
//}
//foreach (var b in VelocityInlets)
//{
//    if (b.setVelState.GetIsNeedToCal(taskIndex))
//    {
//        setUVW_x(b.velocity, b.indexNumber, taskIndex);
//    }
//}
//T getT<T>() where T : class, new()
//{
//    if (velocityInterpolation == null || velocityInterpolation.GetType() != typeof(T))
//        return new T();
//    else
//        return velocityInterpolation as T;
//}
//if (isRhieChow)
//{
//    if (AppRes.computeCenter.isTransient)
//    {
//        TransientRhieChowInterpolation transient = getT<TransientRhieChowInterpolation>();
//        transient.Create(this);
//        transient.SetRhieChow(pressureCorInSIMPLE, pressureInSIMPLE);
//        velocityInterpolation = transient;
//    }
//    else
//    {
//        RhieChowInterpolation rhiechow = getT<RhieChowInterpolation>();
//        rhiechow.Create(this);
//        rhiechow.SetRhieChow(pressureCorInSIMPLE, pressureInSIMPLE);
//        velocityInterpolation = rhiechow;
//    }
//}
//else
//{
//    VelocityInterpolationBase velbase = getT<VelocityInterpolationBase>();
//    velbase.Create(this);
//    velocityInterpolation = velbase;
//}